COVID-19 Data Analysis

ITALY REGIONS OVERVIEW

Click on a region to remove

Double-click on a region to isolate

Data repository: link

Jupyter Notebook repository: link

REGIONS DATA TABLE

In [1]:
import copy
import json
import requests
import datetime as dt

import numpy as np
import pandas as pd

from matplotlib import pyplot as plt
from matplotlib.dates import date2num, num2date
from matplotlib import dates as mdates
from matplotlib import ticker
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

import plotly.graph_objects as go
import plotly.offline as pyo

from scipy.interpolate import interp1d
import scipy.stats as sps

import pyjags
import arviz as az

from IPython.display import clear_output, HTML, Markdown

pyo.init_notebook_mode()
%matplotlib inline

Rt

Italy regions $R_t$ (effective reproduction number in time $t$) calculated with Bayesian approach (Bettencourt & Ribeiro 2008 and Kevin Systrom 2020).

(read here for details)

In [2]:
defined_p = [.95, .5]
In [3]:
# find and fix outliers using Hampel filter
# Impl from: https://towardsdatascience.com/outlier-detection-with-hampel-filter-85ddf523c73d

def hampel_filter_pandas(input_series, window_size, n_sigmas=3.0):

    k = 1.4826 # scale factor for Gaussian distribution
    new_series = input_series.copy()

    # helper lambda function 
    MAD = lambda x: np.median(np.abs(x - np.median(x)))
    
    # the use of min_periods is to have rolling window extend towards
    # the end of the data series; in effect, we can apply hampel filter
    # to most recent observations
    # taken from: https://stackoverflow.com/questions/48953313/pandas-rolling-window-boundary-on-start-end-of-series/48953314#48953314
    rolling_window_size = 2*window_size+1
    rolling_median = input_series.rolling(
        window=rolling_window_size,
        min_periods=(rolling_window_size//2),
        center=True).median()
    rolling_mad = k * input_series.rolling(
        window=rolling_window_size,
        min_periods=(rolling_window_size//2),
        center=True).apply(MAD)
    # print(f'rolling_mad = {rolling_mad}, rolling_median = {rolling_median}')
    diff = np.abs(input_series - rolling_median)
    
    where = diff > (n_sigmas * rolling_mad)
    indices = np.argwhere(where.to_numpy()).flatten()
    new_series[indices] = rolling_median[indices]
    
    return new_series, indices
In [5]:
np.seterr(all='ignore')

states = None
url = 'https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv'
states = pd.read_csv(url,
                     usecols=['data', 'denominazione_regione', 'nuovi_positivi'],
                     parse_dates=['data'],
                     index_col=['denominazione_regione', 'data'],
                     squeeze=True).sort_index()

# We create an array for every possible value of Rt
R_T_MAX = 12
r_t_range = np.linspace(0, R_T_MAX, R_T_MAX*100+1)

# Gamma is 1/serial interval
# https://wwwnc.cdc.gov/eid/article/26/7/20-0282_article
# https://www.nejm.org/doi/full/10.1056/NEJMoa2001316
GAMMA = 1/7.5

def prepare_cases(cases, cutoff=1, filter_win=7, filter_std=2, smooth_win=14, smooth_std=5):

    cases[cases < 0] = 0

    cases, cases_outliers = hampel_filter_pandas(cases, filter_win, filter_std)
    
    smoothed = cases.rolling(smooth_win,
        win_type='gaussian',
        min_periods=1,
        center=True).mean(std=smooth_std).round()
    smoothed[smoothed <= 0] = 0 #1e-10
    
    #idx_start = np.searchsorted(smoothed, cutoff)
    
    smoothed = smoothed.iloc[cutoff:]
    original = cases.loc[smoothed.index]
    
    return original, smoothed

def highest_density_interval(pmf, p=.9):
    # If we pass a DataFrame, just call this recursively on the columns
    if(isinstance(pmf, pd.DataFrame)):
        return pd.DataFrame([highest_density_interval(pmf[col], p=p) for col in pmf],
                            index=pmf.columns)
    
    cumsum = np.cumsum(pmf.values)
    
    # N x N matrix of total probability mass for each low, high
    total_p = cumsum - cumsum[:, None]
    
    # Return all indices with total_p > p
    lows, highs = np.array([]), np.array([])
    j = 0
    while not lows.size or not highs.size:
        lows, highs = (total_p > (p-j)).nonzero()
        j += .05
        if j > 1:
            break
    
    if not len(lows):
        print("no lo")
        print(pmf)
        print(total_p)
        lows = np.array([0])
    if not len(highs):
        print("no hi")
        highs = np.array([1])
    
    try:
        # Find the smallest range (highest density)
        best = (highs - lows).argmin()
    except Exception as e:
        print("ERR {}".format(e))
        best = 0
    
    try:
        low = pmf.index[lows[best]]
    except Exception as err:
        print(f"lo{p} ERROR: {err}")
        low = 0
    try:
        high = pmf.index[highs[best]]
    except Exception as err:
        print(f"hi{p} ERROR: {err}")
        high = 1
    
    return pd.Series([low, high],
                     index=[f'Low_{p*100:.0f}',
                            f'High_{p*100:.0f}'])

def HDI_of_grid(probMassVec, credMass=0.95):
    sortedProbMass = np.sort(probMassVec, axis=None)[::-1]
    HDIheightIdx = np.min(np.where(np.cumsum(sortedProbMass) >= credMass))
    HDIheight = sortedProbMass[HDIheightIdx]
    HDImass = np.sum(probMassVec[probMassVec >= HDIheight])
    idx = np.where(probMassVec >= HDIheight)[0]
    return {'indexes':idx, 'mass':HDImass, 'height':HDIheight}


def HDI_of_grid_from_df(pmf, p):
    # If we pass a DataFrame, just call this recursively on the columns
    if(isinstance(pmf, pd.DataFrame)):
        return pd.DataFrame([HDI_of_grid_from_df(pmf[col], p=p) for col in pmf],
                            index=pmf.columns)
    res = HDI_of_grid(pmf, p)
    #print(res["indexes"])
    lo_idx = res["indexes"][0]
    hi_idx = res["indexes"][-1]
    
    lo = pmf.index[lo_idx]
    hi = pmf.index[hi_idx]
    
    return pd.Series([lo, hi],
                     index=[f'Low_{p*100:.0f}',
                            f'High_{p*100:.0f}'])


def HDIs(pmf, P=[.95, .5]):
    RES = []
    for p in P:
        res = HDI_of_grid_from_df(pmf, p=p)
        RES.append(res)
    return RES

def get_posteriors(sr, sigma=0.15):

    # (1) Calculate Lambda
    lam = sr[:-1].values * np.exp(GAMMA * (r_t_range[:, None] - 1))

    
    # (2) Calculate each day's likelihood
    likelihoods = pd.DataFrame(
        data = sps.poisson.pmf(sr[1:].values, lam),
        index = r_t_range,
        columns = sr.index[1:])
    likelihoods[likelihoods==0] = 1e-15
    
    # (3) Create the Gaussian Matrix
    process_matrix = sps.norm(loc=r_t_range,
                              scale=sigma
                             ).pdf(r_t_range[:, None]) 

    # (3a) Normalize all rows to sum to 1
    process_matrix /= process_matrix.sum(axis=0)
    
    # (4) Calculate the initial prior
    prior0 = sps.gamma(a=4).pdf(r_t_range)
    prior0 /= prior0.sum()

    # Create a DataFrame that will hold our posteriors for each day
    # Insert our prior as the first posterior.
    posteriors = pd.DataFrame(
        index=r_t_range,
        columns=sr.index,
        data={sr.index[0]: prior0}
    )
    
    # We said we'd keep track of the sum of the log of the probability
    # of the data for maximum likelihood calculation.
    log_likelihood = 0.0

    # (5) Iteratively apply Bayes' rule
    for previous_day, current_day in zip(sr.index[:-1], sr.index[1:]):

        #(5a) Calculate the new prior
        current_prior = process_matrix @ posteriors[previous_day]
        
        #(5b) Calculate the numerator of Bayes' Rule: P(k|R_t)P(R_t)
        numerator = likelihoods[current_day] * current_prior
        
        #(5c) Calcluate the denominator of Bayes' Rule P(k)
        denominator = np.sum(numerator)
        
        # Execute full Bayes' Rule
        posteriors[current_day] = numerator/denominator
        
        # Add to the running sum of log likelihoods
        #p_err = np.seterr(all='raise')
        #try:
        #    log_likelihood += np.log(denominator)
        #except Exception as err:
        #    print(f"log err {err} (d={denominator})")
        #    log_likelihood += np.log(1+1e-10)
        #np.seterr(**p_err)
    
    return posteriors, log_likelihood

# targets = ~states.index.get_level_values('denominazione_regione')
# states_to_process = states.loc[states]

#_results = {}

########################### SIGMAS ###########################
sigmas = np.linspace(1/20, 1, 20)

#for state_name, cases in states.groupby(level='denominazione_regione'):
#    print(f"log likelihoods for {state_name}...")
#    smoothed = []
#    cutoff = 20
#    while not len(smoothed):
#        new, smoothed = prepare_cases(cases, cutoff=cutoff)
#        cutoff -= 1
#    
#    result = {}
#    
#    # Holds all posteriors with every given value of sigma
#    result['posteriors'] = []
#    
#    # Holds the log likelihood across all k for each value of sigma
#    result['log_likelihoods'] = []
#    
#    for sigma in sigmas:
#        posteriors, log_likelihood = get_posteriors(smoothed, sigma=sigma)
#        result['posteriors'].append(posteriors)
#        result['log_likelihoods'].append(log_likelihood)
#    
#    # Store all results keyed off of state name
#    _results[state_name] = result
#    #clear_output(wait=True)

########################### SIGMAS ###########################

sigma = .25

results = None

for name, cases in states.groupby(level='denominazione_regione'):
    #print(f"Rt {name}...")
    original, smoothed = prepare_cases(cases)
    #if name == "Sardegna":
    #    print(original)
    #    print(smoothed)
    #    plt.plot(original.values, ls=":", c="k")
    #    plt.plot(smoothed.values, ls="-", c="r")

    #sigmas = np.linspace(1/20, 1, 20)
    
    ## Each index of this array holds the total of the log likelihoods for
    ## the corresponding index of the sigmas array.
    #total_log_likelihoods = np.zeros_like(sigmas)
    #
    ## Loop through each state's results and add the log likelihoods to the running total.
    #for state_name, result in results.items():
    #total_log_likelihoods += _results[name]['log_likelihoods']
    #
    ## Select the index with the largest log likelihood total
    #max_likelihood_index = total_log_likelihoods.argmax()
    #
    ## Select the value that has the highest log likelihood
    #sigma = sigmas[max_likelihood_index]
    #print(f"sigma = {sigma}")
    
    # Note that we're fixing sigma to a value just for the example
    posteriors, log_likelihood = get_posteriors(smoothed, sigma=sigma)

    # Note that this takes a while to execute - it's not the most efficient algorithm
    #hdis = highest_density_interval(posteriors, p=defined_p)
    HDIS = HDIs(posteriors, P=defined_p)

    most_likely = posteriors.idxmax().rename('ML')

    # Look into why you shift -1
    result = most_likely
    for hdis in HDIS:
        result = pd.concat([result, hdis], axis=1)

    results = pd.concat([results, result])

#clear_output()
In [6]:
def plot_rt(result, ax, state_name, P=[.95, .5]):
    
    ax.set_title(f"{state_name}")
    
    # Colors
    ABOVE = [1,0,0]
    MIDDLE = [1,1,1]
    BELOW = [0,.75,0]
    cmap = ListedColormap(np.r_[
        np.linspace(BELOW,MIDDLE,25),
        np.linspace(MIDDLE,ABOVE,25)
    ])
    color_mapped = lambda y: np.clip(y, .5, 1.5)-.5
    
    index = result['ML'].index.get_level_values('data')
    values = result['ML'].values
    update = results.index.get_level_values('data')[-1].date()
    
    # Plot dots and line
    ax.plot(index, values, c='k', zorder=1, alpha=.25)
    ax.scatter(index,
               values,
               s=20,
               lw=.25,
               c=cmap(color_mapped(values)),
               edgecolors='k', zorder=2)
    
    extended = pd.date_range(start=pd.Timestamp('2020-02-24'),
                         end=index[-1]+pd.Timedelta(days=1))
    
    # Aesthetically, extrapolate credible interval by 1 day either side
    for j, p in enumerate(P):
        lowfn = interp1d(date2num(index),
                         result[f'Low_{p*100:.0f}'].values,
                         bounds_error=False,
                         fill_value='extrapolate')

        highfn = interp1d(date2num(index),
                          result[f'High_{p*100:.0f}'].values,
                          bounds_error=False,
                          fill_value='extrapolate')



        ax.fill_between(extended,
                        lowfn(date2num(extended)),
                        highfn(date2num(extended)),
                        color='k',
                        alpha=.1*(j+1),
                        lw=0,
                        zorder=3+j)

    ax.axhline(1.0, c='k', lw=1, label='$R_t=1.0$', alpha=.25);
    
    # Formatting
    ax.xaxis.set_major_locator(mdates.MonthLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%b'))
    ax.xaxis.set_minor_locator(mdates.DayLocator())
    
    ax.yaxis.set_major_locator(ticker.MultipleLocator(1))
    ax.yaxis.set_major_formatter(ticker.StrMethodFormatter("{x:.1f}"))
    ax.yaxis.tick_right()
    ax.spines['left'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.margins(0)
    ax.grid(which='major', axis='y', c='k', alpha=.1, zorder=-2)
    ax.margins(0)
    ax.set_ylim(0.0, 5.0)
    ax.set_xlim(pd.Timestamp('2020-02-24'), result.index.get_level_values('data')[-1]+pd.Timedelta(days=1))
    plt.suptitle(
        f"Italian Regions Real Time $R_t$ up to {update} (HDI {','.join([f'{p:.0%}' for p in sorted(P)])})", 
        fontsize=30, y=1.02)
    fig.set_facecolor('w')
    return values[-1]
In [7]:
_R0 = {}
ncols = 3
nrows = 7

fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(15, nrows*3))

_ = """smr = results.rolling(7,
        win_type='gaussian',
        min_periods=1,
        center=True).mean(std=2)"""

for i, (state_name, result) in enumerate(results.groupby('denominazione_regione')):
    last_Rt = plot_rt(result, axes.flat[i], state_name, P=defined_p)
    _R0.update({state_name: last_Rt})

fig.tight_layout()
fig.set_facecolor('w')

with open("R0_regioni.json", "w") as f:
    json.dump(_R0, f)
In [7]:
results.to_pickle("Rt-regions-time-series.pkl")
In [8]:
BO_COLOR = [.7,.7,.7]
NO_COLOR = [179/255,35/255,14/255,.75]
OK_COLOR = [.5,.5,.5]
KO_COLOR = [35/255,179/255,14/255,.75]
ERROR_BAR_COLOR = [.3,.3,.3]

def plot_standings(msr, figsize=None, p=.9, days=2):
    display(Markdown("***"))
    mr = msr.groupby(level=0)[['ML', f'High_{p*100:.0f}', f'Low_{p*100:.0f}']].tail(days)
    mr = mr.groupby("denominazione_regione")[['ML', f'High_{p*100:.0f}', f'Low_{p*100:.0f}']].mean()
    display(HTML(mr.to_html()))
    mr.sort_values('ML', inplace=True)
    update = results.index.get_level_values('data')[-1].date()
    
    if not figsize:
        figsize = ((15.9/50)*len(mr)+.1,2.5)
        
    fig, ax = plt.subplots(figsize=figsize)

    if days == 1:
        txt = "last day"
    else:
        txt = f"last {days} days mean up to {update}"
    ax.set_title(f'Most Recent $R_t$ by Region ({txt}) - C.I. {p:.0%}')
    
    err = mr[[f'Low_{p*100:.0f}', f'High_{p*100:.0f}']].sub(mr['ML'], axis=0).abs()
    bars = ax.bar(np.arange(mr.index.size),
                  mr['ML'],
                  width=.825,
                  color=BO_COLOR,
                  ecolor=ERROR_BAR_COLOR,
                  capsize=5,
                  error_kw={'alpha':1, 'lw':1},
                  yerr=err.values.T)

    for bar, ml, h90 in zip(bars, mr["ML"], mr[f"High_{p*100:.0f}"]):
        if ml > 1:
            bar.set_color(NO_COLOR)
        #if ml < .75:
        #    bar.set_color(OK_COLOR)
        if ml == 1:
            bar.set_color(BO_COLOR)
        if ml < 1:
            bar.set_color(KO_COLOR)

    labels = mr.index.unique().values  #.replace({'District of Columbia':'DC'})
    ax.set_xticks(np.arange(mr.index.size))
    ax.set_xticklabels(labels, rotation=90, fontsize=11)
    ax.margins(0)
    ax.set_ylim(0,max(mr[f"High_{p*100:.0f}"])*1.1)
    ax.axhline(1.0, linestyle=':', color='k', lw=1)

    ax.legend(handles=[
                        Patch(label='Likely out of control', color=NO_COLOR),
                        #Patch(label='Unsure', color=BO_COLOR),
                        #Patch(label='Near 0', color=OK_COLOR),
                        Patch(label='Likely under control', color=KO_COLOR)
                    ],
                    title='$R_t$ level',
                    ncol=2,
                    loc='upper left',
                    columnspacing=.75,
                    handletextpad=.5,
                    handlelength=1)

    #leg._legend_box.align = "left"
    fig.set_facecolor('w')
    plt.show()
    plt.close(fig="all")
    return fig, ax

plot_standings(results, figsize=(15,5), p=min(defined_p), days=1);
plot_standings(results, figsize=(15,5), p=min(defined_p), days=7);

ML High_50 Low_50
denominazione_regione
Abruzzo 1.21 1.63 0.75
Basilicata 1.08 1.75 0.42
Calabria 0.94 1.42 0.44
Campania 1.41 1.65 1.17
Emilia-Romagna 1.31 1.56 1.06
Friuli Venezia Giulia 1.12 1.53 0.68
Lazio 1.07 1.30 0.84
Liguria 1.30 1.65 0.94
Lombardia 1.27 1.47 1.06
Marche 0.87 1.31 0.42
Molise 0.60 1.46 0.21
P.A. Bolzano 0.92 1.38 0.44
P.A. Trento 1.49 2.05 0.91
Piemonte 1.34 1.64 1.04
Puglia 1.35 1.69 1.01
Sardegna 1.27 1.58 0.95
Sicilia 0.92 1.25 0.57
Toscana 1.31 1.60 1.02
Umbria 1.34 1.76 0.90
Valle d'Aosta 0.64 1.83 0.21
Veneto 1.04 1.28 0.79

ML High_50 Low_50
denominazione_regione
Abruzzo 1.270000 1.714286 0.797143
Basilicata 1.221429 1.938571 0.475714
Calabria 1.034286 1.520000 0.504286
Campania 1.621429 1.872857 1.370000
Emilia-Romagna 1.448571 1.707143 1.188571
Friuli Venezia Giulia 1.270000 1.705714 0.815714
Lazio 1.335714 1.568571 1.101429
Liguria 1.362857 1.728571 0.988571
Lombardia 1.482857 1.691429 1.267143
Marche 0.922857 1.365714 0.452857
Molise 0.591429 1.465714 0.207143
P.A. Bolzano 1.107143 1.587143 0.582857
P.A. Trento 1.400000 1.984286 0.774286
Piemonte 1.382857 1.692857 1.067143
Puglia 1.455714 1.808571 1.101429
Sardegna 1.572857 1.900000 1.244286
Sicilia 1.004286 1.341429 0.654286
Toscana 1.424286 1.725714 1.121429
Umbria 1.494286 1.947143 1.030000
Valle d'Aosta 0.622857 1.942857 0.192857
Veneto 1.202857 1.447143 0.951429

GROWTH

Read here for current method PDF.

In [8]:
states = None
url = 'https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv'
states = pd.read_csv(url,
                     usecols=['data', 'denominazione_regione', 
                              'totale_positivi', 'dimessi_guariti', 'deceduti'
                             ],
                     parse_dates=['data'],
                     index_col=['denominazione_regione', 'data'],
                     squeeze=True).sort_index()

days = states.index.get_level_values('data').unique()

with open("priors.npy", "rb") as f:
    bartlett_samples = np.load(f)

fig, ax = plt.subplots(3, 7, figsize=(15,7))

for i, state in enumerate(states.index.get_level_values('denominazione_regione').unique()):
    f = states.loc[state]
    I = f['totale_positivi']
    R = f['dimessi_guariti']
    D = f['deceduti']
    ax.flat[i].set_title(state, fontsize=10)
    ax.flat[i].fill_between(
        days,
        0, I, color="b"
    )
    ax.flat[i].fill_between(
        days,
        I, R+I, color="lightgreen"
    )
    ax.flat[i].fill_between(
        days,
        R+I, R+I+D, color="r"
    )
    ax.flat[i].set_xticks([])
    ax.flat[i].set_yticks([])
In [9]:
def slope(x, y0, yz):
    alpha = (yz - y0) / (x.size - 1)
    return alpha

def exp_unif(x, beta, y0, yz):
    with np.errstate(all='ignore'):
        y = np.nan_to_num((np.exp(beta * x) - 1)/(np.exp(beta) - 1))
        y = y * (yz - y0) + y0
        
    alpha = slope(x, y0, yz)

    return (y, alpha)

window = 14
x = np.linspace(0, 1, window)
In [10]:
fig, ax = plt.subplots(figsize=(5,2))
ax.text(0, .1, r"$\times$", color="k", fontsize=20, ha='left', va='center')
ax.text(0, .3, r"$\to$", color="b", fontsize=20, ha='left', va='center')
ax.text(0, .5, r"$\downarrow$", color="g", fontsize=20, ha='left', va='center')
ax.text(0, .7, r"$\nearrow$", color="y", fontsize=20, ha='left', va='center')
ax.text(0, .9, r"$\uparrow$", color="r", fontsize=20, ha='left', va='center')

ax.text(.1, .1, "Data ERROR", color="k", fontsize=20, ha='left', va='center')
ax.text(.1, .3, "None (constant)", color="b", fontsize=20, ha='left', va='center')
ax.text(.1, .5, "Exp. Decay", color="g", fontsize=20, ha='left', va='center')
ax.text(.1, .7, "Almost Linear", color="y", fontsize=20, ha='left', va='center')
ax.text(.1, .9, "Exponential", color="r", fontsize=20, ha='left', va='center')

ax.set_xticks([])
ax.set_yticks([])

ax.set_title("Symbols Legenda")
plt.show()
plt.close(fig="all")
In [11]:
beta_lim = 487

modelString = f"""
model {{
  for ( r in 1:R ) {{
    # inverse gamma variance
    tau[r] ~ dgamma( 0.001 , 0.001 )
    sigma[r] <- 1 / sqrt( tau[r] )
    # beta priors
    beta[r] ~ dunif( -{beta_lim} , {beta_lim} )
    
    for ( t in 1:W ) {{
      # define lim(beta->0)
      H[r,t] <- ( exp( beta[r] * x[t] ) - 1 )
      f[r,t] <- ifelse( nu[r] == 1 , x[t] , H[r,t] / nu[r] )
      # expected and observed
      E[r,t] <- f[r,t] * (yw[r] - y1[r]) + y1[r]
      y[r,t] ~ dnorm( E[r,t] , tau[r] )
    }}
    
    # define lim(beta->0)
    nu[r] <- ifelse( exp( beta[r] ) - 1 == 0 , 1 , exp( beta[r] ) - 1 )
  }}
}}
"""

df = pd.read_csv(
    "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv",
     usecols=['data', 'denominazione_regione', 'totale_casi', 'deceduti'],
     parse_dates=['data'],
     index_col=['denominazione_regione', 'data'],
     squeeze=True).sort_index()

countries = df.index.get_level_values("denominazione_regione").unique().sort_values()
days = df.index.get_level_values("data").unique()

w_min = -window-1
w_max = -1

growth = {}

for what in df.columns:
    display(Markdown(f"## {what}"))
    growth.update({
        (what, 'x'): {},
        (what, 'growth'): {},
        (what, 'beta'): {},
        (what, 'alpha'): {},
    })
    y = np.array([])
    Y = np.array([])
    beta_lims = np.array([])
    for country in countries:
        growth[(what, 'x')].update({country: None})
        growth[(what, 'growth')].update({country: None})
        growth[(what, 'beta')].update({country: "-"})
        growth[(what, 'alpha')].update({country: "-"})
        T = df.loc[country][what][w_min:w_max].values
        _T = copy.deepcopy(T)
        if _T[0] > _T[-1]:
            _T = _T[0] * np.ones(window)
        if not y.size:
            Y = copy.deepcopy(T)
            y = copy.deepcopy(_T)
        else:
            Y = np.vstack([Y, T])
            y = np.vstack([y, _T])
        
    model_data = {
        "R": countries.size,
        "W": window,
        "x": x,
        "y": y,
        "y1": y[:,0],
        "yw": y[:,window-1],
    }
    jags_posterior_model = pyjags.Model(
        code=modelString, 
        data=model_data, 
        chains=4, 
        adapt=2000,
        progress_bar=False
    )

    #_ = jags_posterior_model.sample(2000, vars=["beta", "sigma"])
    
    jags_posterior_samples = jags_posterior_model.sample(8000, vars=["beta", "sigma"])
    
    pyjags_data = az.from_pyjags(jags_posterior_samples)

    summary = az.summary(pyjags_data, var_names="^beta", filter_vars="regex", hdi_prob=.99, kind="stats")
    summary_diag = az.summary(pyjags_data, var_names="^beta", filter_vars="regex", hdi_prob=.99, kind="diagnostics")
    
    sigma_smry = az.summary(pyjags_data, var_names="^sigma", filter_vars="regex", hdi_prob=.99, kind="stats")
    sigma_smry_diag = az.summary(pyjags_data, var_names="^sigma", filter_vars="regex", hdi_prob=.99, kind="diagnostics")
    
    #display(Markdown(f"### sigma"))
    
    _ = """for i, par in enumerate(sigma_smry.index):
        axs = az.plot_trace(
            pyjags_data.posterior.sigma[:,:,i],
            figsize=(20,2)
        )
        for j, ax in enumerate(axs):
            ax[j].set_title(f"sigma\n{countries[i]}")
        plt.show()
        plt.close(fig="all")
    
    axs = az.plot_density(
        pyjags_data.posterior.sigma,
        hdi_prob=.995,
        point_estimate="mode",
        figsize=(20,2*7)
    )
    for i, ax in enumerate(axs):
        ax.set_title(f"sigma\n{countries[i]}")
    plt.show()
    plt.close(fig="all")"""
    
    _ = """ax = az.plot_forest(
        pyjags_data, combined=True,
        var_names="^sigma",
        filter_vars="regex",
        hdi_prob=.995,
        figsize=(12,6),
    )
    ax[0].axvline(0, c="k", alpha=.2)
    ax[0].set_yticklabels(countries[::-1])
    ax[0].set_title(fr"Estimated $\sigma$ for '{what}' (99.5% HPD)")
    plt.show()
    plt.close(fig="all")
    
    display(HTML(sigma_smry.set_index(countries).to_html()))"""
    #display(Latex(sigma_smry_diag.set_index(countries).to_latex(longtable=True, )))
    
    display(Markdown(f"### beta"))
    
    _ = """for i, par in enumerate(summary.index):
        axs = az.plot_trace(
            pyjags_data.posterior.beta[:,:,i],
            figsize=(20,2)
        )
        for j, ax in enumerate(axs):
            ax[j].set_title(f"beta\n{countries[i]}")
        plt.show()
        plt.close(fig="all")
    
    axs = az.plot_density(
        pyjags_data.posterior.beta,
        hdi_prob=.995,
        figsize=(20,2*7)
    )
    for i, ax in enumerate(axs):
        ax.set_title(f"beta\n{countries[i]}")
    plt.show()
    plt.close(fig="all")"""

    beta_smry = summary[(summary["hdi_0.5%"]>-beta_lim*.5) & (summary["hdi_99.5%"]<beta_lim*.5)]
    
    ax = az.plot_forest(
        pyjags_data, combined=True,
        var_names="^beta",
        filter_vars="regex",
        hdi_prob=.995,
        figsize=(12,6),
    )
    ax[0].set_xlim(beta_smry["hdi_0.5%"].min(), beta_smry["hdi_99.5%"].max())
    ax[0].axvline(0, c="k", alpha=.2)
    ax[0].set_yticklabels(countries[::-1])
    ax[0].set_title(fr"Estimated $\beta$ parameter growth for '{what}' (99.5% HPD)")
    plt.show()
    plt.close(fig="all")
    
    display(HTML(summary.set_index(countries).to_html()))
    #display(Latex(summary_diag.set_index(countries).to_latex(longtable=True, )))
    
    display(Markdown(f"### fitted"))
    fig, ax = plt.subplots(7, 3, figsize=(15, 3*7), sharex=True)
    for par in range(countries.size):
        y1 = Y[par,0]
        yw = Y[par,-1]
        
        mu = summary.loc[f"beta[{par}]"]["mean"]
        lo = summary.loc[f"beta[{par}]"]["hdi_0.5%"]
        hi = summary.loc[f"beta[{par}]"]["hdi_99.5%"]
        sd = summary.loc[f"beta[{par}]"]["sd"]
        
        prmu, alpha = exp_unif(x, mu, y1, yw)
        prlo, _ = exp_unif(x, lo, y1, yw)
        prhi, _ = exp_unif(x, hi, y1, yw)
        
        bartlett_s, bartlett_p = sps.bartlett(
            bartlett_samples, 
            pyjags_data.posterior.beta[0,:,par].values.flatten(),
            pyjags_data.posterior.beta[1,:,par].values.flatten(),
            pyjags_data.posterior.beta[2,:,par].values.flatten(),
            pyjags_data.posterior.beta[3,:,par].values.flatten()
        )
        
        if yw < y1:
            growth[(what, 'x')][countries[par]] = -1
            growth[(what, 'growth')][countries[par]] = "ERR"
            txt = r"$\times$"
            col = "k"
        elif bartlett_p >= .05:
            # constant
            growth[(what, 'x')][countries[par]] = 0
            growth[(what, 'growth')][countries[par]] = "None"
            growth[(what, 'beta')][countries[par]] = 0
            growth[(what, 'alpha')][countries[par]] = 0
            txt = r"$\to$"
            col = "b"
        else:
            growth[(what, 'beta')][countries[par]] = f"{mu:.2f}"
            growth[(what, 'alpha')][countries[par]] = f"{alpha:.2f}"
            if lo <= 0 <= hi:
                # linear
                growth[(what, 'x')][countries[par]] = 2
                growth[(what, 'growth')][countries[par]] = "lin"
                txt = r"$\nearrow$"
                col = "y"
            # exp or dec
            elif mu < 0:
                # dec
                growth[(what, 'x')][countries[par]] = 1
                growth[(what, 'growth')][countries[par]] = "dec"
                txt = r"$\downarrow$"
                col = "g"
            else:
                # exp
                growth[(what, 'x')][countries[par]] = 3
                growth[(what, 'growth')][countries[par]] = "exp"
                txt = r"$\uparrow$"
                col = "r"
        
        ax.flat[par].plot(days[w_min:w_max], Y[par], 'x:', c="k", ms=5, label="observed")
        ax.flat[par].plot(days[w_min:w_max], prmu)
        ax.flat[par].fill_between(
            days[w_min:w_max], prlo, prhi, alpha=.2, label="fitted"
        )
        ax.flat[par].text(
            .9, .1, txt, color=col,
            fontsize=20, ha='center', va='center', transform=ax.flat[par].transAxes,
            bbox={'facecolor': 'w', 'alpha': .75, 'edgecolor': 'k'}
        )
        ax.flat[par].plot(
            days[w_min:w_max], alpha * np.arange(window) + y1, c="k", alpha=.2
        )
        ax.flat[par].set_title(fr"{countries[par]} $\bar{{\beta}}={mu}$ $\alpha={alpha:.2f}$")
        ax.flat[par].xaxis.set_major_locator(mdates.WeekdayLocator())
        ax.flat[par].xaxis.set_major_formatter(mdates.DateFormatter('%d %b'))
        ax.flat[par].xaxis.set_minor_locator(mdates.DayLocator())
        ax.flat[par].legend(fontsize=8, loc="upper left")
    plt.show()
    plt.close(fig="all")

growth_df = pd.DataFrame(growth)

display(Markdown(f"## Summary"))
fig, ax = plt.subplots(figsize=(15,2))

ax.plot(growth_df[('deceduti', 'x')],
    "+-", lw=0, ms=14,
    label="deaths"
)
ax.plot(growth_df[('totale_casi', 'x')],
    "x-", lw=0, ms=12,
    label="cases"
)
ax.set_title(f"Latest estimated growths (last {window} days)")
ax.set_xlim(-1, countries.size+3)
ax.set_ylim(-1.5, 3.5)
ax.set_xticks(np.arange(0, countries.size))
ax.set_xticklabels(countries, rotation=90)
ax.set_yticks([-1,0,1,2,3])
ax.set_yticklabels(["ERR", "None", "log", "lin", "exp"])
ax.fill_between(
    np.arange(-1, countries.size+1), -1.5, -0.5, color="k", alpha=.1
)
ax.fill_between(
    np.arange(-1, countries.size+1), -0.5, 0.5, color="b", alpha=.1
)
ax.fill_between(
    np.arange(-1, countries.size+1), 0.5, 1.5, color="g", alpha=.1
)
ax.fill_between(
    np.arange(-1, countries.size+1), 1.5, 2.5, color="y", alpha=.1
)
ax.fill_between(
    np.arange(-1, countries.size+1), 2.5, 3.5, color="r", alpha=.1
)
ax.legend()
plt.show()
plt.close(fig="all");

display(Markdown('Growths summary:'))

display(HTML(growth_df.to_html(columns=[
    ('deceduti','growth'), ('deceduti','beta'), ('deceduti','alpha'),
    ('totale_casi','growth'), ('totale_casi','beta'), ('totale_casi','alpha'),
])))

smry = pd.DataFrame(
    {
        'ERR': {
            'deaths': growth_df[('deceduti', 'x')][growth_df[('deceduti', 'x')]==-1].count(),
            'cases': growth_df[('totale_casi', 'x')][growth_df[('totale_casi', 'x')]==-1].count(),
         },
        'None': {
            'deaths': growth_df[('deceduti', 'x')][growth_df[('deceduti', 'x')]==0].count(),
            'cases': growth_df[('totale_casi', 'x')][growth_df[('totale_casi', 'x')]==0].count(),
         },
        'dec': {
            'deaths': growth_df[('deceduti', 'x')][growth_df[('deceduti', 'x')]==1].count(),
            'cases': growth_df[('totale_casi', 'x')][growth_df[('totale_casi', 'x')]==1].count(),
         },
        'lin': {
            'deaths': growth_df[('deceduti', 'x')][growth_df[('deceduti', 'x')]==2].count(),
            'cases': growth_df[('totale_casi', 'x')][growth_df[('totale_casi', 'x')]==2].count(),
         },
        'exp': {
            'deaths': growth_df[('deceduti', 'x')][growth_df[('deceduti', 'x')]==3].count(),
            'cases': growth_df[('totale_casi', 'x')][growth_df[('totale_casi', 'x')]==3].count(),
         },
    }
)

display(HTML(smry.to_html()))

display(Markdown(fr"""
For cumulative total cases

- {smry.loc["cases"]["exp"]} regions {smry.loc["cases"]["exp"]/21:.2%} are likely exponentially growing
- {smry.loc["cases"]["lin"]} regions {smry.loc["cases"]["lin"]/21:.2%} are likely almost linearly growing
- {smry.loc["cases"]["dec"]} regions {smry.loc["cases"]["dec"]/21:.2%} are likely exponentially decaying
- {smry.loc["cases"]["None"]} regions {smry.loc["cases"]["None"]/21:.2%} are likely constant

For cumulative total deaths:

- {smry.loc["deaths"]["exp"]} regions {smry.loc["deaths"]["exp"]/21:.2%} are likely exponentially growing
- {smry.loc["deaths"]["lin"]} regions {smry.loc["deaths"]["lin"]/21:.2%} are likely almost linearly growing
- {smry.loc["deaths"]["dec"]} regions {smry.loc["deaths"]["dec"]/21:.2%} are likely exponentially decaying
- {smry.loc["deaths"]["None"]} regions {smry.loc["deaths"]["None"]/21:.2%} are likely constant


In {smry.loc["cases"]["exp"]/21:.0%} regions, cases
are growing exponetially and in {(smry.loc["cases"]["dec"]+smry.loc["cases"]["None"])/21:.0%}
the situation is likely solved or about to resolution (no growth or exponential decay).

On the other hand, in {(smry.loc["deaths"]["dec"]+smry.loc["deaths"]["None"])/21:.0%} regions
deaths are not growing or decaying and {smry.loc["deaths"]["exp"]/21:.0%}
are showing an exponential growth of deaths.

Linearly growing regions ({smry.loc["cases"]["lin"]/21:.0%} for cases and
{smry.loc["deaths"]["lin"]/21:.0%} for deaths) situation is not severe but likely
not currently resolving.
"""))

deceduti

beta

mean sd hdi_0.5% hdi_99.5%
denominazione_regione
Abruzzo -0.569 282.230 -486.846 476.544
Basilicata 1.926 280.684 -479.958 484.038
Calabria 1.685 281.324 -486.885 477.562
Campania -2.147 0.490 -3.620 -0.876
Emilia-Romagna 8.600 2.765 3.945 18.110
Friuli Venezia Giulia 0.492 0.819 -1.928 2.874
Lazio 0.416 0.254 -0.263 1.154
Liguria -2.103 1.083 -5.694 0.646
Lombardia -1.374 0.324 -2.316 -0.499
Marche -0.300 281.888 -483.293 481.016
Molise 0.613 280.736 -477.642 485.211
P.A. Bolzano 0.356 280.395 -477.362 486.469
P.A. Trento -0.410 281.252 -479.173 484.658
Piemonte -1.939 0.344 -2.959 -1.036
Puglia -2.151 1.410 -6.103 0.767
Sardegna -0.473 280.567 -477.236 486.876
Sicilia 1.021 11.006 -45.857 74.693
Toscana -0.290 1.161 -3.730 2.770
Umbria -1.931 280.857 -486.333 476.933
Valle d'Aosta 1.535 282.170 -477.727 486.629
Veneto 0.631 0.363 -0.382 1.645

fitted

totale_casi

beta

mean sd hdi_0.5% hdi_99.5%
denominazione_regione
Abruzzo -0.914 0.188 -1.427 -0.399
Basilicata -1.289 0.532 -2.812 0.156
Calabria 1.849 0.214 1.250 2.441
Campania 2.095 0.099 1.827 2.372
Emilia-Romagna 0.061 0.085 -0.171 0.299
Friuli Venezia Giulia 0.935 0.154 0.514 1.359
Lazio 1.581 0.050 1.441 1.718
Liguria 1.934 0.243 1.294 2.614
Lombardia -0.132 0.090 -0.378 0.113
Marche 0.830 0.168 0.371 1.308
Molise 2.264 0.531 0.909 3.853
P.A. Bolzano -0.298 0.225 -0.926 0.305
P.A. Trento -1.244 0.162 -1.683 -0.797
Piemonte 0.604 0.086 0.364 0.838
Puglia 0.645 0.122 0.317 0.984
Sardegna 0.796 0.222 0.169 1.396
Sicilia 0.857 0.130 0.489 1.218
Toscana 0.921 0.111 0.619 1.231
Umbria 0.868 0.136 0.488 1.251
Valle d'Aosta -2.610 0.752 -4.983 -0.757
Veneto 0.644 0.139 0.261 1.030

fitted

Summary

Growths summary:

deceduti totale_casi
growth beta alpha growth beta alpha
Abruzzo None 0 0 dec -0.91 12.46
Basilicata None 0 0 lin -1.29 2.46
Calabria None 0 0 exp 1.85 6.00
Campania dec -2.15 0.31 exp 2.10 19.85
Emilia-Romagna exp 8.60 12.46 lin 0.06 48.31
Friuli Venezia Giulia lin 0.49 0.15 exp 0.94 7.38
Lazio lin 0.42 0.46 exp 1.58 33.08
Liguria lin -2.10 0.08 exp 1.93 16.92
Lombardia dec -1.37 1.69 lin -0.13 80.23
Marche None 0 0 exp 0.83 10.62
Molise None 0 0 exp 2.26 0.77
P.A. Bolzano None 0 0 lin -0.30 5.23
P.A. Trento None 0 0 dec -1.24 1.69
Piemonte dec -1.94 0.62 exp 0.60 28.69
Puglia lin -2.15 0.08 exp 0.65 16.15
Sardegna None 0 0 exp 0.80 5.92
Sicilia lin 1.02 0.15 exp 0.86 35.23
Toscana lin -0.29 0.15 exp 0.92 25.00
Umbria None 0 0 exp 0.87 5.62
Valle d'Aosta None 0 0 dec -2.61 0.62
Veneto lin 0.63 1.54 exp 0.64 75.15
ERR None dec lin exp
deaths 0 10 3 7 1
cases 0 0 3 4 14

For cumulative total cases

  • 14 regions 66.67% are likely exponentially growing
  • 4 regions 19.05% are likely almost linearly growing
  • 3 regions 14.29% are likely exponentially decaying
  • 0 regions 0.00% are likely constant

For cumulative total deaths:

  • 1 regions 4.76% are likely exponentially growing
  • 7 regions 33.33% are likely almost linearly growing
  • 3 regions 14.29% are likely exponentially decaying
  • 10 regions 47.62% are likely constant

In 67% regions, cases are growing exponetially and in 14% the situation is likely solved or about to resolution (no growth or exponential decay).

On the other hand, in 62% regions deaths are not growing or decaying and 5% are showing an exponential growth of deaths.

Linearly growing regions (19% for cases and 33% for deaths) situation is not severe but likely not currently resolving.


In [12]:
json_regions = "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-json/dpc-covid19-ita-regioni.json"
with requests.get(json_regions) as req:
    data = json.loads(req.content.decode('utf-8-sig'))
In [13]:
print("FIRST ENTRY DATE: {}".format(
    data[0]["data"]
    )
)
print("LAST  ENTRY DATE: {}".format(
    data[-1]["data"]
    )
)
period = (
    dt.datetime.strptime(data[-1]["data"], "%Y-%m-%dT%H:%M:%S") -
    dt.datetime.strptime(data[0]["data"], "%Y-%m-%dT%H:%M:%S")
).days

print("COVERAGE: {} days".format(period))
print("CURRENT DATE IS: {}".format(dt.datetime.now().strftime("%Y-%m-%d %H:%M:%S")))
FIRST ENTRY DATE: 2020-02-24T18:00:00
LAST  ENTRY DATE: 2020-08-17T17:00:00
COVERAGE: 174 days
CURRENT DATE IS: 2020-08-17 23:34:11
In [14]:
x = []   # datetime x array
# _x = []  # integer x array
# yC = []  # new confirmed cases array
# yD = []  # new deaths array
# yR = []  # new recovered array
# yP = []  # new infected array

TOTyC = {}  # confirmed cases array
TOTyD = {}  # deaths array
TOTyR = {}  # recovered array
TOTyP = {}  # infected array
TOTyr = {}  # mortality rate
TOTyk = {}  # recovery rate

TOTyPric = {}  # ricoverati
TOTyPint = {}  # intensiva
TOTyPiso = {}  # isolamento
In [15]:
for entry in data:

    # x values
    date = dt.datetime.strptime(entry["data"], "%Y-%m-%dT%H:%M:%S")
    if date not in x:
        x.append(date)
        
    region = entry["denominazione_regione"]
    if region not in TOTyC:
        TOTyC.update({region: np.array([])})
        TOTyD.update({region: np.array([])})
        TOTyR.update({region: np.array([])})
        TOTyP.update({region: np.array([])})
        TOTyr.update({region: np.array([])})
        TOTyk.update({region: np.array([])})
        TOTyPric.update({region: np.array([])})
        TOTyPint.update({region: np.array([])})
        TOTyPiso.update({region: np.array([])})
    # y TOT values
    TOTyC[region] = np.append(TOTyC[region], entry["totale_casi"])
    TOTyD[region] = np.append(TOTyD[region], entry["deceduti"])
    TOTyR[region] = np.append(TOTyR[region], int(entry["dimessi_guariti"]))
    TOTyP[region] = np.append(TOTyP[region], entry["totale_positivi"])
    TOTyPric[region] = np.append(TOTyPric[region], entry["ricoverati_con_sintomi"])
    TOTyPint[region] = np.append(TOTyPint[region], entry["terapia_intensiva"])
    TOTyPiso[region] = np.append(TOTyPiso[region], entry["isolamento_domiciliare"])
    
    if entry["totale_casi"]:
        TOTyr[region] = np.append(TOTyr[region], entry["deceduti"] / entry["totale_casi"])
        TOTyk[region] = np.append(TOTyk[region],
                                  (int(entry["dimessi_guariti"])) / 
                                  entry["totale_casi"])
    else:
        TOTyr[region] = np.append(TOTyr[region], .0)
        TOTyk[region] = np.append(TOTyk[region], .0)

_C = [TOTyC[r][0] for r in TOTyC]
minC = max(_C)
_D = [TOTyD[r][0] for r in TOTyD]
minD = max(_D)
_R = [TOTyR[r][0] for r in TOTyR]
minR = max(_R)
_P = [TOTyP[r][0] for r in TOTyP]
minP = max(_P)

_ric = [TOTyPric[r][0] for r in TOTyPric]
minric = min(_ric)
_int = [TOTyPint[r][0] for r in TOTyPint]
minint = min(_int)
_iso = [TOTyPiso[r][0] for r in TOTyPiso]
miniso = min(_iso)

TOTAL CASES

In [16]:
fig = go.Figure()

for region in TOTyC:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyC[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (cases)", "xanchor": "center", "x": 0.5},
    yaxis_title="cases",
)

pyo.iplot(fig)
In [17]:
fig = go.Figure()

for region in TOTyC:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyC[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (cases, log scale)", "xanchor": "center", "x": 0.5},
    yaxis_title="cases", yaxis_type="log",
)

pyo.iplot(fig)

DEATHS

In [18]:
fig = go.Figure()

for region in TOTyD:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyD[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (deaths)", "xanchor": "center", "x": 0.5},
    yaxis_title="deaths",
)

pyo.iplot(fig)
In [19]:
fig = go.Figure()

for region in TOTyD:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyD[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (deaths, log scale)", "xanchor": "center", "x": 0.5},
    yaxis_title="deaths", yaxis_type="log",
)

pyo.iplot(fig)

RECOVERED

In [20]:
fig = go.Figure()

for region in TOTyR:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyR[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (recovered)", "xanchor": "center", "x": 0.5},
    yaxis_title="recovered",
)

pyo.iplot(fig)
In [21]:
fig = go.Figure()

for region in TOTyR:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyR[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (recovered, log scale)", "xanchor": "center", "x": 0.5},
    yaxis_title="recovered",  yaxis_type="log",
)

pyo.iplot(fig)

INFECTED

In [22]:
fig = go.Figure()

for region in TOTyP:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyP[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (infected)", "xanchor": "center", "x": 0.5},
    yaxis_title="infected",
)

pyo.iplot(fig)
In [23]:
fig = go.Figure()

for region in TOTyP:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyP[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (infected, log scale)", "xanchor": "center", "x": 0.5},
    yaxis_title="infected",  yaxis_type="log",
)

pyo.iplot(fig)

MORTALITY RATE

!!! PLEASE NOTE !!!

These rates are only useful for SIRD epidemiological model (read here for details) not to define COVID-19 actual rates.

In [24]:
fig = go.Figure()

for region in TOTyr:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyr[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="cross",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.2%}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696', "tickformat": ',.0%',},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (mortality rate)", "xanchor": "center", "x": 0.5},
    yaxis_title="percentage",
)

pyo.iplot(fig)

RECOVERY RATE

In [25]:
fig = go.Figure()

for region in TOTyk:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyk[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="cross",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.2%}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696', "tickformat": ',.0%',},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (recovery rate)", "xanchor": "center", "x": 0.5},
    yaxis_title="percentage",
)

pyo.iplot(fig)

OSPITALIZED DETAILS

In [26]:
fig = go.Figure()

for region in TOTyPric:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyPric[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (symptomatic hospitalized)", "xanchor": "center", "x": 0.5},
    yaxis_title="number",
)

pyo.iplot(fig)
In [27]:
fig = go.Figure()

for region in TOTyPint:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyPint[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (in emergency)", "xanchor": "center", "x": 0.5},
    yaxis_title="number",
)

pyo.iplot(fig)
In [28]:
fig = go.Figure()

for region in TOTyPiso:
    fig.add_trace(go.Scatter(
        x=x, y=TOTyPiso[region],
        mode='lines+markers',
        marker_size=3, marker_symbol="circle",
        line_shape='spline',
        name=region,
        hovertemplate="%{text}"+
        "<br>%{x}"+
        "<br>%{y:.0f}",
        text=[region for _ in range(len(x))]
    ))

fig.update_layout(legend_orientation="h",
    showlegend=True, plot_bgcolor='rgba(0,0,0,0)', 
    yaxis={"gridcolor": '#bdbdbd', "zerolinecolor": '#969696'},
    xaxis={"gridcolor": '#bdbdbd'},
    title={"text": "ITALY regions (home isolated)", "xanchor": "center", "x": 0.5},
    yaxis_title="number",
)

pyo.iplot(fig)
In [ ]:
 

© 2020 Max Pierini & Sandra Mazzoli & Alessio Pamovio

Exported from italy/regions-overview.ipynb committed by Max Pierini on Mon Aug 31 10:59:12 2020 revision 1, 8eb7ef8